homework 10, version 0

Submission by: Jazzy Doe (jazz@mit.edu)

Homework 10: Climate modeling II

18.S191, fall 2020

Let's create a package environment:

In Lecture 23 (video above), we looked at a 2D ocean model that included two physical processes: advection (flow of heat) and diffusion (spreading of heat). This homework includes the model from the lecture, and you will be able to experiment with it yourself!

The model is written in a way that it can be extended with more physical processes. In this homework we will add two more effects, introduced in the Energy Balance Model from our last homework: absorbed and emitted radiation.

Exercise 1 - Advection-diffusion

Included below is the two-dimensional advection-diffusion model from Lecture 23. To keep this homework concise, we have only included the code. To see the original notebook with more detailed comments, use the link below:

Advection & diffusion

Notice that both functions have a main method with the following signature:

(::Array{Float64,2}, ::ClimateModel) maps to ::Array{Float64,2}.

As we will see later, ClimateModel contains the grid, the velocity vector field and the simulation parameters.

Data structures

Let's look at our first type, Grid. Notice that it only has one 'constructor function', which takes N (number of longitudinal grid points) and L (longitudinal size in meters) as arguments. The struct contains more fields, these are precomputed and stored for performance.

Next, let's look at three types.

Two structs: OceanModel and OceanModelParameters, and an abstract type: ClimateModel.

Timestepping

The OceanModel struct contains a complete description of the model being simulated. The next struct, ClimateModelSimulation, contains a model, together with the simulation results. It is mutable: timestepping the model means modifying a ClimateModelSimulation object.

Exercise 1.1 - Running the model

In the next few cells, we set up a simulation. We have included an interactive visualisation of the simulation.

👉 Familiarize yourself with the simulation through interaction. Get a sense for each parameter by changing their values.

Uncomment (Ctrl+/ or Cmd+/) one of the lines below to choose between the different velocity fields:

Choose the initial temperature state:

We define our ocean simulation. Run this cell again to reset the simulation to the initial state.

Simulation controls

Click to show the velocity field or to show temperature anomalies instead of absolute values

Vary the current speed U: 1.0 [× reference]

Vary the diffusivity κ:

10000.0

[m²/s]

Velocity field for a single circular vortex

Quasi-realistic ocean velocity field $\overrightarrow{u}=(u,v)$

Our velocity field is given by an analytical solution to the classic wind-driven gyre problem, which is given by solving the fourth-order partial differential equation:

$$-{ϵ}_M{\hat{\mathrm{\nabla }}}^4\hat{\mathrm{\Psi }}+\frac{\partial \hat{\mathrm{\Psi }}}{\partial \hat{x}}=\mathrm{\nabla }\times \hat{\tau }\mathbf{z},$$

where the hats denote that all of the variables have been non-dimensionalized and all of their constant coefficients have been bundles into the single parameter ${ϵ}_M\equiv {\displaystyle \frac{\nu }{\beta L^3}}$.

The solution makes use of an advanced asymptotic method (valid in the limit that $ϵ\ll 1$) known as boundary layer analysis (see MIT course 18.305 to learn more).

Some simple initial temperature fields

👉 Some parameters have a physical meaning (κ, u and v), while other parameters control our numerical process. Choose two of these numerical parameters, and describe their effect on the simulation's runtime and the simulation's accuracy.

Exercise 2 - Complexity

In this class we have restricted ourself to small problems ($N_t<100$ timesteps and $N_{x;y}<30$ spatial grid-cells) so that they can be run interactively on an average computer. In state-of-the-art climate modelling however, the goal is to push the numerical resolution $N$ to be as large as possible (the grid spacing $\mathrm{\Delta }t$ or $\mathrm{\Delta }x$ as small as possible), to resolve physical processes that improve the realism of the simulation (see below).

Here, we provide a simple estimate of the computational complexity of climate models, which reveals a substantial challenge to the improvement of climate models.

Our climate model algorithm can be summarized by the recursive formula:

$$T_{i,j}^{n+1}=T_{i,j}^n+\mathrm{\Delta }t\ast \left(\text{tendencies}\right)$$

for each time step $n\in \{1,\dots ,M\}$, and for each grid point $i\in \{1,\dots ,N_x\},j\in \{1,\dots ,N_y\}$.

Our goal is to simulate an ocean of fixed size (e.g. $6000$ km), for a fixed amount of time (e.g. $100$ years). By choosing smaller $\mathrm{\Delta }t$, $\mathrm{\Delta }x$ and $\mathrm{\Delta }y$, we get a more accurate simulation, but for the same size, it requires more steps. i.e.

$$M=\mathcal{O}(\mathrm{\Delta }{t}^{-1})N_x=\mathcal{O}(\mathrm{\Delta }{x}^{-1})N_y=\mathcal{O}(\mathrm{\Delta }{y}^{-1}).$$

Now, the total runtime of our simulation is proportional to the number of steps we need to take, which is $M\cdot N_x\cdot N_y$. For a fixed aspect ratio $N_y=2N_x$, we get

$$\text{runtime}=\mathcal{O}(M)\mathcal{O}(2N_x^2)=\mathcal{O}(M)\mathcal{O}(N_x^2).$$

Exercise 2.1

For constant $M$, we want to verify that $\text{runtime}=\mathcal{O}(N_x^2)$ holds for our numerical model.

👉 Write a function model_runtime that takes N as an argument, and sets up a model with grid of resolution N, and returns the runtime of a single timestep!.

👉 Call your runtime function on a range of values for N, and use a plot to demonstrate that the predicted runtime complexity holds.

Exercise 2.2 - The CFL condition on $\mathrm{\Delta }t$

In Exercise 1, look for the definition of Δt. It is currently set to 12*60*60 (12 hours).

👉 Double Δt and run the simulation again. You should see that it runs faster, great! Now, keep doubling Δt until you see something 'strange'. Describe what you see.

What you experienced is a numerical instability of the discretization method in our simulation. This is not caused by floating point errors – it is a theoretical limitation of our method.

To ensure the stability of our finite-difference approximation for advection, heat should not be displaced more than one grid cell in a single timestep. Mathematically, we can ensure this by checking that the distance $L_{CFL}\equiv max(|\overrightarrow{u}|)\mathrm{\Delta }t$ is significantly less than the width $\mathrm{\Delta }x=\mathrm{\Delta }y$ of a single grid cell:

$$L_{CFL}\equiv max(|\overrightarrow{u}|)\mathrm{\Delta }t\ll \mathrm{\Delta }x$$

or

$$\mathrm{\Delta }t\ll \frac{\mathrm{\Delta }x}{max(|\overrightarrow{u}|)},$$

which is known as the Courant-Freidrichs-Levy (CFL) condition.

The exact meaning of $\ll$ here depends on the simulation at hand, but in our case, it means that there is at least an order of magnitude difference:

$$\mathrm{\Delta }t<0.1\frac{\mathrm{\Delta }x}{max(|\overrightarrow{u}|)}.$$

Given below is a function CFL_advection that takes a ClimateModel and computes the CFL value: $\mathrm{\Delta }x/max(|\overrightarrow{u}|)$.

👉 Using the interactive simulation of Exercise 1, verify that the CFL condition is (somewhat) true. Increase the magnitude of $\overrightarrow{u}$ until you start to see numerical artefacts. Now, halve the value fo $\mathrm{\Delta }t$, and again, increase the magnitude of $\overrightarrow{u}$. You should find that halving $\mathrm{\Delta }t$ allows for twice the velocity magnitude. The same link exists between $\mathrm{\Delta }t$ and $\mathrm{\Delta }x$.

The CFL inequality states that if we want to decrease the grid spacing $\mathrm{\Delta }x$ (or increase the resolution $N_x$), we also have to decrease the timestep $\mathrm{\Delta }t$ by the same factor. In other words, the timestep can not be thought of as fixed – it depends on the spatial resolution: $\mathrm{\Delta }t\equiv \mathrm{\Delta }{t}_0\mathrm{\Delta }x$. This means that $M=\mathcal{O}(N_x)$.

Revisiting our complexity equation, we now have

$$\text{runtime}=\mathcal{O}(M)\mathcal{O}(N_x^2)=\mathcal{O}(N_x^3).$$

In other words, in a 2-D model, 2x the spatial resolution requires 8x the computational power.

Exercise 2.3 - Moore's Law

In practice, state-of-the-art climate models are 3-D, not 2-D. It turns out that to preserve the aspect ratio of oceanic motions, the vertical grid resolution should also be increased $N_z\propto N_x$, such that in reality the computational complexity of climate models is:

$$\text{runtime}=\mathcal{O}(N_x^4).$$

This is the fundamental challenge of high-performance climate computing: to increase the resolution of the models by a factor of $2$, the model's run-time increases by a factor of $2^4=16$.

The figure below shows how the grid spacing of state-of-the-art climate models has decreased from $500$ km in 1990 (FAR) to $100$ km in the 2010s (AR4). In other words, grid resolution increased by a factor of $5$ in 20 years.

Moore's law is the observation that the number of transistors in a dense integrated circuit doubles about every two years. In the context of climate modelling, we can interpret this as meaning that the computational complexity allowed by our best high-performance computers $\mathcal{C}$ doubles every two years:

$$\mathcal{C}(t)=\mathcal{C}(2020)\ast 2^{(t-2020)/2}$$

Present-day simulations have a grid spacing of $\mathrm{\Delta }x=30$ km (or about $N_x=L/\mathrm{\Delta }x\approx \frac{20000\text{ km}}{30\text{ km}}\approx 700$).

👉 By extrapolating Moore's law forward into the future, estimate how long it would take for $\mathrm{\Delta }x$ to reach to the $500$ meter scale of clouds, one of the important climate processes ($N_x=L/\mathrm{\Delta }x\approx 40000$).

Exercise 3 - Radiation

In Homework 9, we used a zero-dimensional (0-D) Energy Balance Model (EBM) to understand how Earth's average radiative imbalance results in temperature changes:

$$\begin{array}{rlr}C\frac{\partial T}{\partial t}=& \phantom{+}\frac{S(1-\alpha )}{4}& \text{(absorbed radiation)}\\ & -(A-BT)& \text{(outgoing radiation)}\end{array}$$

This week we will do the same, but now in two-dimensions (2-D), where in addition to heat being added or removed from the system by radiation, heat can be transported around the system by oceanic advection and diffusion (see also Lectures 22 & 23 for 1-D and 2-D advection-diffusion). The governing equation for temperature $T(x,y,t)$ in our coupled climate model is:

$$\begin{array}{rlr}\frac{\partial T}{\partial t}=& \phantom{+}u(x,y)\frac{\partial T}{\partial x}+v(x,y)\frac{\partial T}{\partial y}& \text{(advection)}\\ & +\kappa (\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2})& \text{(diffusion)}\\ & +\frac{S(x,y)(1-\alpha (T))}{4C}& \text{(absorbed radiation)}\\ & -\frac{(A-BT)}{C}.& \text{(outgoing radiation)}\end{array}$$

Parameters

Below we define two new types: RadiationOceanModel and RadiationOceanModelParameters. Notice the similarities and differences with our original model. By making RadiationOceanModel also a subtype of ClimateModel, we will be able to re-use much of our original code.

Notice that this struct has Base.@kwdef in front of its definition. This allows us to:

1. assign default values to the struct fields, directly inside the definition, and

2. automatically create an easier constructor function that takes keyword arguments:

Exercise 3.1 - Absorbed radiation

The ocean absorbs solar radiation, increasing the temperature. Just like in our EBM model, we model the ocean as a surface with a temperature-dependent albedo: $\alpha (T)$. A high albedo means that more sunlight is reflected, and less is absorbed. We model the albedo function as:

An area of ocean below 0°C is covered in ice, which is more reflective, and therefore absorbs less solar radiation. In our EBM model, this positive feedback leads to a bifurcation: under the same external conditions, the climate system has multiple equilibria.

In this week's two-dimensional model, the factor $\alpha$ is also two-dimensional: instead of a global albedo, every grid cell has its own temperature, which determines its own albedo, $\alpha (T(x,y,t))$. We can now have an ocean with warm and cold regions, which absorb different amounts of radiation. The same positive feedback can have a local effect.

Here is a second method for α that takes a 2D array T with the current ocean temperatures and a RadiationOceanModel, and returns the 2D array of albedos. We use the dot operator to apply α pointwise to T, also called broadcasting.

Solar insolation

Our rectangular grid represents the North Atlantic Ocean, stretching from the equator (latitude = 0°) to the North Pole (latitude = 90°) in the latitudinal direction (y), and from the east coast of North America to the west coast of Europe in the longitudinal direction (x). In reality, climate models have to explicitly deal with the curvature of the Earth when constructing their model grids. Here, we will just treat the North Atlantic Ocean as a rectangle with roughly the correct dimensions.

Just like the albedo, every grid cell will have a local amount of solar insolation. In our model, we use the annual average at the latitude of a grid cell $S(y)/4$. This is given by: (source)

Note that latitude is the horizontal axis in this graph, not the vertical.

👉 Write a method absorbed_solar_radiation that takes a 2D array T with the current ocean temperatures and a RadiationOceanModel, and returns the tendencies corresponding to absorbed radiation. This is the analogue of advect and diffuse.

Exercise 3.2 - Outgoing radiation

Just like in our EBM from before, when our ocean heats up by absorbing solar radiation, it also emits radiation back to space. This outgoing thermal radiation is what allows the ocean to eventually come to an equilibrium. The difference here is that each individual grid cell of our model radiatives according to it's local temperature $T_{i,j}$.

👉 Write a method outgoing_termal_radiation that takes a 2D array T with the current ocean temperatures and a RadiationOceanModel, and returns the tendencies corresponding to outgoing radiation. This is the analogue of advect and diffuse.

Exercise 3.3 - Running the model

Let's define a new timestep! method for our new model. This is exactly the same as our advection-diffusion model, with the addition of the two radiation tendencies.

We can now simulate our radiation ocean model, reusing much of the code from our advection-diffusion simulation.

👉 Play around with the simulation to find the effect of each parameter. In particular, discover the effect of the solar insolation, S_mean, the initial temperatures, T_init, the liquid and ice albedos: α0 and αi, and the amount of emitted heat at 0°C, A.

Exercise 3.4 - Stable states

So far, we are able to set up a model and run it interactively. You see that the model quickly goes from the initial temperatures to a stable state: a state with balanced energy (radiation out, radiation in). Changing the initial state slightly will probably result in the same stable state.

But let's see what happens when we initialize with extremely high or low initial temperatures...

👉 For $S=1380$ (present-day value) and default parameters, does T_init_value=-50 give a different result than T_init_value=+50? What about T_init_value=+55? By trying various values for T_init_value, how many stable states do you find?

👉 Answer the same question for $S=1000$ and $S=2000$.

Exercise 3.5

That's right, we have found a bifurcation! Under the right conditions, there are multiple stable states. But under different conditions, there is only one stable state.

Hypothesis: The multiple stable states at $S=1380$ are caused by the feedback effect of the ice albedo. A cold ocean reflects more heat and stays cold, a warm ocean absorbs more heat and stays warm.

👉 Find a way to confirm this hypothesis by changing the model parameters in the interactive model. Describe your findings.

👉 Why is the number of stable states different for a lower or higher value of $S$?

Exercise 4 - Bifurcation diagram

So far, we are able to set up a model and run it interactively, and we discovered that by changing the initial value of S_mean, we find a different number of stable states.

In this final exercise, we will generate a visualization to help us understand the relationship between S_mean and the number of stable states. Instead of running a single model interactively, we write a function that takes the model parameters as input, runs the model until equilibrium, and returns the final mean temperature. We will run this high-level function for various initial values, generating a single graph: the bifurcation diagram.

Exercise 4.1 - Equilibrium temperature

👉 Write a function eq_T that takes two argments, S and T_init_value, that sets up a radiation ocean model with S as S_mean, and with T_init_value as the constant initial temperature. Run the model until you have reached equilibrium (approximately), and return the average temperature.

👉 Use the function eq_T to verify your answer to Exercise 3.4.

Exercise 4.2

👉 Make a bifurcation diagram. For a couple of pairs (S,T_init), calculate the equilibrium temperature. Create a scatter plot with $S$ on the horizontal axis, and $T_{\text{equilibrium}}$ on the vertical.

This calculation can take quite some time. Some tips:

1. Start out with a small amount of points.

2. Run the eq_T calculations in one cell, and store the result as a vector. Generate the scatter plot in a different cell.

3. You can replace a map with ThreadsX.map to run it in parallel across multiple CPU cores. More on this below.

Parallelize a map

You use map to apply a function to each element of a vector. Why not use a for loop? One nice property of map code is that you only describe the operation for each element, not the order to run the operations in. This means that it can be easily parallelized by the computer. For example, on a computer with 4 cores, the computation map(sqrt, 1:100) can be parallelized by handling 1:25 on the first core, 26:50 on the second, etc., at the same time.

In Julia, many functional primitives (map, filter, sum, maximum, all, and more) have an automatically multithreaded version, in the ThreadsX.jl package. As a demo, compare the two cells below. (You can see the runtime in the bottom right of a cell.) You might need to run the cells a second time, the first time includes Julia's compiler doing its thing.

Exercise XX: Lecture transcript

(MIT students only)

Please see the link for the transcript document on Canvas. We want each of you to correct about 500 lines, but don’t spend more than 20 minutes on it. See the the beginning of the document for more instructions.

👉 Please mention the name of the video(s) and the line ranges you edited:

Appendix

Function library

Just some helper functions used in the notebook.